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Abstract. Two-fluid ideal plasma equations are a generalized form of the ideal MHD equa- 
tions in which electrons and ions are considered as separate species. The design of efficient 
numerical schemes for the these equations is complicated on account of their non-linear na- 
ture and the presence of stiff source terms, especially for high charge to mass ratios and for 
low Larmor radii. In this article, we design entropy stable finite difference schemes for the 
two-fluid equations by combining entropy conservative fluxes and suitable numerical diffusion 
operators. Furthermore, to overcome the time step restrictions imposed by the stiff source 
terms, we devise time-stepping routines based on implicit-explicit (IMEX)-Runge Kutta (RK) 
schemes. The special structure of the two-fluid plasma equations is exploited by us to design 
IMEX schemes in which only local (in each cell) linear equations need to be solved at each 
time step. Benchmark numerical experiments are presented to illustrate the robustness and 
accuracy of these schemes. 



1. Introduction 

An ensemble of plasma consists of ions, electrons and neutral particles. These particles interact 
through both short range (e.g. collisions) and long range ( e.g. electromagnetic) forces. Plasmas 
are increasingly used in spacecraft propulsion, controlled nuclear fusion and in circuit breakers 
in the electrical power industry. 

Under the assumption of quasi-neutrality ( i.e. charge density difference between ions and 
electrons is neglected), the flow of plasmas is modeled by the ideal MHD equations (see [8]). 
Although, the ideal MHD equations have been successfully employed in modeling and simulating 
plasma flows, this model is derived by ignoring the Hall effect and treating plasma flows as 
single fluid flows. These effects are very important for many applications, e.g. space plasmas, 
Hall current thrusters, field reversal configurations for magnetic plasma confinement and for fast 
magnetic reconnection. 

In this article, we consider the more general ideal two- fluid model (see [H] , [E] , [2] ) for collision- 
less plasmas. In the ideal two-fluid equations, electrons and ions are treated as different fluids by 
allowing them to posses different velocities and temperatures. Assuming local thermodynamical 
equilibrium, we write the two- fluid equations as a system of balance laws (see [§]): 

(1.1) <9 t u + divf(u) =s(u), (x, ()6l 3 x (0, oo). 

Here, u = u(x, y, z, t) is the vector of non-dimensional conservative variables, 

(1-2) u = {pi,piVi,Ei,p e ,p e v e ,E e ,B,E, (j), ip} T . 

Here, the subscripts {i, e} refer to the ion and electron species respectively, P{i, e } are t ne densities, 
v {i.e} = { v \i e}i v \i e }> v {i e}) are ^ ne velocities, -E{j, e } are the total energies, B = (B x , B V ,B Z ) is 
the magnetic field, E = (E x , E v , E z ) is the electric field and (f>,ip are the potentials. The flux 
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vector f = (f x , f v , P) can be written as, 

fi(Uj) 



(1.3) f(u) = 

with a £ {i, e}, and 



fe(Ue) 
f*m (^m ) 



where f Q (u a ) 
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for any vector K = {K x , K v , K z ). Here u Q = {p a , p a v a , E a } T , a E {i, e}, u m = {B, E, <j>, tp} T , 
P{i,e} are t ne pressures, £, k are penalizing speeds (see |14j ) and c = c/uf 1 is the normalized speed 
of light. Also, vf is the reference thermal velocity of ion. Writing the flux in component form 
(see (|1.3p . (|1.4l0 . we observe that the first two components of the flux, f Q (u Q ), a £ {«,e}, are 
the nonlinear ion and electron Euler fluxes and the third component is the linear Maxwell flux. 
So, the homogeneous part of (jl.lj) is hyperbolic. 
The source term s is given by, 



iftCE + Vi x B) 
9 £pi(E.Vi) 


-4p-,9 e (E + v e x B) 

( L5 ) S(U)=< ! -^ e (E.V e ) 









h r e p e v e ) 

TePe) 



with the charge to mass ratios r Q = q a /m a , a 6 {«,e} and the ion-electron mass ratio A m = 
mi/m e . Two physically significant parameters appear in the source term namely, the normalized 

T 

Larmor radius f — — q .^ XQ and the ion Debye length (normalized with respect to the 

Here, Bo is the reference magnetic field, Eq is 



Larmor radius) X d = A d /r g = ^ e vf 2 /n qi/r g 
the permittivity of free space and xq is the reference length. The ion mass rrii and ion charge 
qi are assumed to be 1. In addition, we assume that both the ion and the electron satisfies the 
ideal gas law: 

Pa 1 



(1.6) 



7 



1 



a e {i, e}, 



with gas constant 7 = 5/3. In the above equations, we use the perfectly hyperbolic form of the 
Maxwell equation (see [11]), which represent the evolution of magnetic field B and electric field 
E. 

The design of numerical schemes for systems of balance laws has undergone rapid devel- 
opment in the last two decades, see [H] for a detailed description of efficient schemes. The 
standard paradigm involves the use of finite volume (conservative finite difference) schemes in 
which the solution is evolved in terms of (approximate) solutions of Riemann problems at cell 
interfaces. Higher order accuracy in space is obtained by non-oscillatory interpolation proce- 
dures of the TVD, ENO and WENO types. An alternative is to use the Discontinuous Galerkin 
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(DG) approach. High-order temporal accuracy results from strong stability preserving (SSP) 
Runge-Kutta (RK) methods. Source terms are included by using operator splitting approaches. 

Although the two-fluid equations are a system of balance laws, standard discretization tech- 
niques may fail to provide a robust approximation. Two major difficulties are present in the 
numerical analysis of the two-fluid equations: I ) the design of suitable numerical fluxes and 2) 
treatment of the source term that becomes increasingly stiffer as more realistic charge to mass 
ratios or more realistic Larmor radii (Debye lengths) are considered. 

Given the above challenges, very few robust numerical schemes exist for the two-fluid equa- 
tions. In |15) . the authors derive a Roe- type Riemann solver. Time updates are performed by 
treating the stiff source term implicitly and the flux terms explicitly. The resulting non-linear 
equations are solved using Newton iterations. This method might be diffusive and may require 
many iterations for each time step. In [5], the authors propose a wave propagation method (see 
[T2"] ) for the spatial discretization. For time updates, a second -order operator splitting approach 
is used. A similar approach is taken in [131 E] , where spatial discretization is based on discon- 
tinuous Galerkin (DG) methods and time update is based on SSP-RK methods. Both of these 
approaches are easy to implement but can be computationally expensive, especially for realistic 
charge to mass ratios. 

Given the state of the art, we propose numerical schemes for the two-fluid equations with the 
following novel features: 

• First, we design entropy stable finite difference discretizations of the two fluid equations. 
The basis of our design is to ensure entropy stability as the two fluid equations satisfy an 
entropy inequality at the continuous level. We use the approach of |17j by constructing 
entropy conservative fluxes and suitable numerical diffusion operators to ensure entropy 
stability. Second-order entropy stable schemes are constructed following the framework 
of a recent paper [6]. 

• We discretize the source term in the two-fluid equations by an IMEX approach: the flux 
terms are discretized explicitly whereas the source term is discretized implicitly. The 
main feature of our schemes is their ability to use the special structure of the two-fluid 
equations that allows us to design IMEX schemes requiring the solution of only local 
(in each cell) linear equations at every time step. This is in contrast to the schemes 
proposed in [15] that required the solution of non-linear iterations. The local equations 
that result from our approach can be solved exactly making our schemes computationally 
inexpensive. 

The rest of this article is organized as follows: In the following Section[2l we obtain an entropy 
estimate for the ideal two- fluid eqns. This result at the continuous level is then used to 

design an entropy stable finite difference scheme in Section [3J In Section 2J we present IMEX- 
RK schemes for the temporal discretization. The resulting, algebraic system of equations is then 
solved exactly. In Section[Sl we simulate the nonlinear soliton propagation in the two-fluid plasma 
and a stiff Riemann problem to demonstrate the robustness and efficiency of these schemes. 

2. Analysis of Continuous Problem 

It is well known that solutions of (jl.ip consists of discontinuities, even for smooth initial data. 
Hence, we need to consider the solutions of in the weak sense. However, uniqueness of the 
solutions is still not guaranteed and we need to supplement (jl.Ij) with additional admissibility 
criteria to obtain a physically meaningful solution. This gives rise to concept of entropy. The 
standard thermodynamic entropies for ion and electron Euler flows are, 

(2.1) e a = ^-^ with s a = log(p a ) - j\og(p a ), ae{i,e}. 
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For the electromagnetic part we consider the quadratic entropy i.e electromagnetic energy, 

r99 , , s |B| 2 + 2 |E| 2 + ^ 2 
(2-2) e m (u m ) = + — . 

We obtain the following entropy estimate, 

Theorem 2.1. Let u = {pi, pi~Vi, E{, p e , p e v e , E e , B, E, <fi, ip} T be a weak solution of the two-fluid 
Eqns. (|l.ip on R 3 x (0, oo). Furthermore, assume that there exist constants p™ ln ,p™ ax and p™ m 
such that, 

< P T n < Pa < PT X , Pa>P™ n >0, «G{«,e}, 

then 

(2.3) f (e t + e e + e m ) dx dy dz <C\ [ (e t + e e + e m ) dx dy dz + C 2 , 

dt J%3 J #3 

wii/i constant C\ and C 2 depending only on p™ m ,/0™ ax , and p™ m . 

Proof. Let us first consider the fluid part of the equations. The entropy fluxes corresponding to 
the flow entropies (12.11) are, 

(2.4) q Q = — = v a e„, a<E{i,e\. 

7-1 

Assuming that u is a smooth solution of , the densities p a and the pressures p a , satisfy, 

dtp a + v Q • Vp a = 0, 

dtPa + 7P«V • V Q + V a • Vp a = 0. 

Using the expression for s a , we get 

d t s a + v Q • Vs Q = 0. 
Combining this with density advection we get entropy conservation, i.e. 

(2.5) d t e a + V • q Q = 0. 

Observe that (|2.5[) implies that the source term does not effect the evolution of fluid entropies. 
For weak solutions, (|2.5[) reduces to entropy inequality, 

(2.6) d t e a + V • q a < 0. 
Integrating over M 3 and adding, 

d f 

(2.7) — / (a + e e )dx dy dz < 0. 

at J M 3 

For controlling the electromagnetic energy, we use the following inequality, 



(2.8) 



/ (pa + \PaV a \ 2 + -Eq) dx dy dz < C 3 / e a dx dy dz + C4, 

JS. 3 JR 3 



for some constants C, C.The proof of (|2.8I) is a simple consequence of the positivity of density 
and pressure and the use of the relative entropy method of Dafermos [5] . We multiply (jl.ip with 
the vector, 

„ E J V1 T 
Oio,B,-,</),-| 

and note that flux terms are still in divergence form. Integrating over the whole space and using 
Cauchy's inequality on the right hand side, we get, 

(2.9) — / e m dx dy dz < C5 [ f e m dx dy dz + f (e^ + e e )dx dy dz \ + Cq. 

Combining it with (|3.22p we obtain (|2.3[) . □ 
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Remark 2.2. Note that above proof of the theorem also gives a bound on the fluid energy of the 
system. 

3. Semi-Discrete Schemes 

In the last section, we showed that solutions of the two-fluid equations satisfy the entropy 
estimate (|2.3p . In this section, we will design (semi-discrete) numerical schemes for the two- fluid 
equations that satisfy a discrete version of the entropy estimate. 

For simplicity, we consider two- fluid eqns. (Il.lj) in two dimensions, i.e., 

(3.1) d t u + dj x (u) + d y {y(u) = s (u). 

We discretize the two dimensional rectangular domain D = (x a ,Xb) x (y a ,yb) uniformly with 
mesh size (Ax, Ay). We define X4 = x a + iAx and yj = y a + jAy, < i < N x , < j < N y , 
such that xjj — x a + N x Ax and yb = y a + N y Ax. The domain is then divided into cells 
Iij = [^-1/2,^+1/2] x [yj-i/2,Vj+i/2\ with x l+1/2 = Xi+Xi+1 and y j+ i/ 2 = Vi+ % j+1 . A standard 
semi-discrete finite difference scheme for the eqn. (|3.ip can be written as, 

(^) + £ (Fr +1/2)j - F tl/2)J -) + ^ (n, +1/2 - F ? ,_ 1/2 ) = S iA U). 

Here, F x +1 ^ 2 ^ and F^ +] y 2 are the numerical fluxes consistent with P and f y respectively, 
and Si.j(U) = s(Ui.j). We introduce the entropy variables V and entropy potential \ k which 
corresponds to the entropy e = {e^, e e , e m } T 

[ v, ) r a Ui e ( (u 4 ) ) 

(3.3) V = <^ V e \ = I 9 Uc e e (u e ) \ , 

where qj^ is the entropy flux for the Maxwell part corresponding to the entropy e m and k 6 {x, y}. 
We will follow the framework of Tadmor ( see [17l [18]) for designing an entropy stable scheme 
for the two-fluid equations. The first step is to design an entropy conservative flux. 

3.1. Entropy conservative flux. We require the following notation: 

[ a ]i+l/2,j = a i+l,j ~ a i,ji a i+l/2,j 
[<Ai,j+l/2 = a i,j+l — a i,ji a i,j + l/2 

Following [17], an entropy conservative flux F = {F X ,F V } is defined as a consistent flux that 
satisfies 

(3-4) [Vf+i/2 /?+ 1/2)j - {x x Wi/2, v [V]^ +1/2 F^. +1/2 = { X y k J+1/2 . 

In general, the relation for conservative flux, Q3.4p provides one equation for several unknowns. 
Hence, entropy conservative numerical flux is not unique. We will now describe entropy conser- 
vative numerical fluxes for the fluid part of the two-fluid equations. 

In [10] . Ismail and Roe have derived an expression for entropy conservative numerical fluxes 
for Euler equations of gas dynamics. As the fluid part of (jl.ip consists of two independent Eulcr 
fluxes, we can use the expression derived in |10) for the entropy conservative numerical flux of 



X k ={ X h e } = { Vjf e fe -q* } 
Y k V T f fe -q fe 

V Am ) V mm Mm ) 



6 



HARISH KUMAR AND SIDDHARTHA MISHRA 



the Euler flows of ion and electron. We need to introduce parametric vectors z a , a £ {i, e}, 



(3.5) 



Pa 



a e {i, e}. 



Hx,2 



Then the entropy conservative numerical flux in x-direction is given by ~F x a i+1 j 2 j = \F a ' i+1/2 j^a i+1/2 j> 



a, i+1/2, j' a,»+l/2,j' a,i+l/2,j 

(3.6) 



] T , with, 



o,i+i/y 



2 5 
z a,i+l/2j" z ai+l/2, : ,-i 



^,2 

■ a,i+l/2,j 

■ a,i+l/2,j 
a,i+l/2,j 



' a,t+l/2,j 



a,i+l/2,j 



- m 



a,j+l/2J r a,t+l/2,j' 



3 _ ^ 

m a,i+l/2,j F a ,i+l/2,j' 



jlX,l 

a,i4 



m a,!+l/2,3 F a,!+l/2j' 



2Z 1 



a, i+1/2 j 



7+1 r a,»+l/2,3 — A X ,2 
„_i [In ^ n.'+l/2,j r ((l i|l/2,j 

, ' Z a, i+1/2, j 



+ z a,i+l/2* a)i+ i/ 2 ,j + Z a.i+l/aj* a ,i+l/2,jJ ' 



and entropy conservative numerical flux in y-dircction is, F y i J+1/2 = [F^] J+1/2 , F^j , j+1/2 , 



ju/,2 



F y 



^^, 1+1/2 ] T , with, 



a,»,j+l/2' a,i,j+l/2' a,i,j+l/2 

(3.7) 



a,i,j+l/2 



— 3 5 

z a,i,j+l/2 z a,i,j+l/2i 



fy,2 

a,i,j+l/2 

■py,3 

a,i,j+l/2 
r a,i,j+l/2 



x a,i,j+l/2 



5 ln 

a, 2 



" t a,iJ+l/2 r a,»,j+l/2' 

m a,iJ+l/2 + m a,i,j+l/2^a,i,j+l/2' 

m 4 TP 4 '' 1 
" t a,iJ+l/2 r a,»,j+l/2' 



1 



7 + 1 f a,ij+l/2 72 

7 _ 1 z l In ^ ^ a,i,j+l/2' a,i j + 1/2 



2Z 1 



a,iJ+l/2 



M/,2 



"a>i,j+l/2 



+ z a,i,j+l/2* ai i,j +1 /2 + z a,i,j+l/2* a> i t j + l/2) ' 

Here, a- +1 / 2 ^ and a[" +1 / 2 denotes the logarithmic means defined as, 

In _ [ a ]»+l/2,j l n _ [ Q kj+l/2 

U ^ ~ [log (a)} i+1/ 2j ' JJ+1/2 " [log (a)] id+1/2 ' 

and 



m. 



_ z a, i+1/2, j 
a,i+l/2,j - =f : 
z a,i+l/2,j 



r __ Z a,i,j+l/2 

l a,i,j+l/2 ~ =f • 
z a,i,j+l/2 



for re {2,3,4,5}. 



Now we will consider the electromagnetic part. Note the Maxwell flux is linear. Then, it is easy 
to check that the entropy conservative numerical flux for the electromagnetic part is 

(3-8) K,i+i,2,i = |(f x (U ro ,ij) + P(U m , 4+1 ,,)), F^. . +1/2 = ^(f(U miiJ ) + P(U ro , iJ+1 )). 
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Combining all the parts, the entropy conservative numerical flux for the Eqn. are given by, 

! F M+l/2j 1 j F L + l/2 

F e,i+l/2j } ' F i,i+l/2 = { F e,i,j+l/2 

F m,i+l/2,j J [ F m,i, j+1/2 

3.2. Numerical diffusion operator. As entropy is dissipated at shocks, we need to add entropy 
stable numerical diffusion operators to avoid spurious oscillations at shocks. Following [18], the 
resulting numerical fluxes arc of the form, 

(3.10) Ff +1/2;j = F* +1/2 - 2 D i+V2[ V ]j+i/2,j' F L+i/2 = F m/2 ~ ^L'+i/al^ta'+i/a- 
with diffusion matrices are given by, 

fo iii Vi x - R x A 1 R xT n v — 7? y A a R yT 

^O.llj iJ J :+l/2 — rL i+l/2j ii i+l/2j rl i+l/2j' 1J ij + l/2 — rt i,j + l/2 ii i,j + l/2 rt i,j + l/2- 

Here R^ x ^ are the matrices of right eigenvectors of Jacobians d u H x,y ^ and A^ x ' y ^ are diagonal 
matrices of eigenvalues in the x- and y-directions, respectively. We will use a Rusanov type 
diffusion operator given by a 18 x 18 matrix, 

= A = diag{(max |A?|)7 5x5 ,( max |Af|)7 5x5 ,( max |Af|)I 8x8 }. 

l<i<5 6<i<10 l<i<18 

We use the eigenvector scaling due to Barth [¥] for defining the eigenvector matrices. 

3.3. Second Order Dissipation Operator. The diffusion operators described above are of 
first order, as the jump term [V] is of order Ax. To obtain the second order accurate scheme, 
we can perform piecewise linear reconstructions of the entropy variable V. However, a straight- 
forward reconstruction of the entropy variables may not be entropy stable. In [BJ, the authors 
have constructed entropy stable second order (and even higher-order) diffusion operators. For 
simplicity, we will consider the diffusion operator, D? +1 / 2 j\Y]i+i/2,j m x-direction only. The 
diffusion operator in y-direction, D^ j+1 ^ 2 [V] ij+1 / 2 can be constructed analogously. We need to 
introduce scaled entropy variables, 

W^ ± = Eg" 1/2ji V iJ . 

If W*'^ are the reconstructed values of W J± in the x-direction, then the corresponding recon- 
structed values P* for V^- are given by, 

K = I"; : w;;-. 

The resulting second order entropy stable flux is then given by, 

( 3 - 12 ) F i+l/2j = F i+l/2 ~ 2" D *+ 1 /2t P;E ] i + 1 /2j' 

where the jump term [P a: ]i+i/2j is given by, 

;p"-;,.:..2.., p;.,.., p?+. 

A sufficient condition for the scheme to be entropy stable (see [BJ) is the existence of a diagonal 
matrix B x > 0, such that, 

[Wli + i/2j=B? +1/2J [W*] i+1/2J , 

i.e. the reconstruction of *W X has to satisfy a sign preserving property along the interfaces of 
each cell. Component-wise this can be written as, 

(3.13) sign([^]) = signal), 

for each component w and w l of and W 1 , respectively. 
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3.4. Reconstruction Procedure. We suppress the j-dependence below for notational conve- 
nience. The reconstruction for 'W x is performed component- wise, so that (|3.13[) is satisfied. Let 
us define jump of component w of the variable ~W X , 

(3.14) k+x/2 = Mm/2- 

Consider <p, a slope limiter satisfying cf)(9^ 1 ) = <f>(9)9~ 1 and define divided differences, 



0r = ^±V£ and fl+ = 6 -fUl. 

Oj-1/2 Oj+1/2 

Then the reconstructed values of w in the cell Ii are 

™i = w i - 2^(^)^-1/2' w+ = w,+ + -^(^ 1 )^ + i/ 2 - 
Using these values we obtain 

[™]i+l/2 = (l - + W+l))) <W' 

This shows that the sign property is satisfied iff 

< 1, V 9 e M. 



In this article, we will use the minmod limiter based reconstruction which satisfies the sign 
preserving property (see [B]). The minmod limiter is given by, 

f 0, if 9 < 0, 

(3.15) c/)(9) = i 0, if < < 1, 

ll, else. 

3.5. Discrete entropy stability. In this section, we prove that scheme given by the flux (13.121) 
is entropy stable i.e. it satisfies a discrete version of the entropy estimate (|2.3[) . We have the 
following result, 

Theorem 3.1. The semi-discrete finite difference scheme (|3.2p , with entropy stable numerical 
flux (|3.12p , is second order accurate for smooth solutions. Furthermore, it satisfies, 

(3.16) y^X e i,i,j + e ^,j + e m ,i,i)AxAy < C 7 y^ y ( e w + + e m ^ 3 )AxAy + C 8 
if conditions for Theorem \2.1\ are satisfied. 

Proof. It is easy to see that the scheme is of second order accuracy, as both the entropy conserva- 
tive flux F and the numerical diffusion operator, are second order accurate for smooth solutions. 
Now, consider the fluid part of (I3.2[) . i.e. 

( 3 ' 1? ) dt 3 + ~Ax ( F «.i+V2,j ~ K,i-l/2,j) + ^ (''"■'.<..; • I 2 ~ F l,i,j ! j) = S a,i,j(U), 

for a € {i, e} with entropy numerical fluxes, 

(3.18) Qf+l/2,j = V i+l/2,j F m/2j _ Xi+l/2,j, Qfj+l/2 = V »J + l/2 F i'j + l/2 ~ Xi,j+l/2- 

Multiplying Q3.17P with Vj^^ and imitating the proof of Theorem 2.2 from [17], we get 

" 2^ ([ V ]7 +1 /2 J D- + 1 /2,,[P a; ] J+ l/2,, + [V]7-l/ 2)i Dt 1/2> ,[Pl i -l/2 > i) 
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^ (jVU i+1/2 Df J+1/2 [P*] iJ+1/2 + [V^.^D^^P*]^^ 

~Ax (^ +1 /2j ~ Qi-l/2,j) _ ^ (Q'L + 1/2 ~ + V I,i,J S «.i.i( U ) 

4^ ([v]7 +1/2J Df +1/2j [p^] l+1/2j + mu^pu^h-i^ 

4^ ([ v ]i + i/ 2 DL- + i/ 2 [ p? 'k J+ i/2 + [V]^- 1/2 Dr, j _ 1/2 [P I/ ]ij-i/ 2 



Here 



Qf+l/2,j - V i+l/2,i F f+l/2,j - Xi+l/2,], and Qf J + 1 /2 - V ij + l/2 F f J + l/2 _ Xij + 1/2- 

are entropy fluxes corresponding to the entropy conservative fluxes F x and F y respectively. Let 
us consider the diffusion terms. Ignoring all the indices, each diffusion term satisfies, 

[V] T D[P] = [V] T i?Ai? T [P] 

= [V] T i?Ai? T (i? T ) ( - 1} [W] 

= (i? T [V]) T A J B([W]) 

= (i? T [V]) T AB(i? T V) 

> 0, 

as both B and A are non-negative diagonal matrices. So, we get 

' 3 + ^ (Qa,»+l/2j ~ QS,i-l/2,j) + ^7 (Qo,iJ+l/2 _ Qa.ij-l/a) - ^Z,i,j S a,i,j • 

A simple calculation shows that, 

V T S ■ ■ — 

v a,i,j a a,y ~~ u - 

This results in the fluid entropy inequality, 

(3 ' 19) d -^T + ^ x ( Q -i + i/2, - QVw)+5j K« + i/ a - ^ °> e 

summing over all the cells we get, 

d 
dt 



Ay 

•1 

(3.20) ^53e aii)J -AaAy<0, ae{i,e}, 

»,i 

Repeating the entropy argument of Dafermos [5] used in Theorem l2.1l we get an discrete energy 
estimate for fluid part, 

(3-21) J2 + l/°*.*J v «Ail 2 + <. (J ) AxAy < C 9 £ e^AxAy + C w , 

i,j id 

Imitating the proof of Theorem 12.11 where integration is replaced by summation, we get, 
(3.22) ^2 e m ,ij Aa;Ay < C n ^](e m ,ij + e Mj + e e> ij)AxAy + Ci 2 . 

Combining with (ET2"0"|) . we get (|3~TB| . □ 

Remark 3.1. In theorem 13. 11 the discrete energy estimate (|3.16p is satisfied only if the electron 
and ion density and pressure (as required by theorem I2.1[) are positive. We assume that this 
positivity holds for the scheme. Currently, it is not possible to prove that this positivity is also a 
consequence of the numerical scheme. However, we expect that the use of positivity preserving 
limiters (like those in |19| ) will enable us to prove positivity. 
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order 




CHI 




Pa 


2 


1 




1 






1/2 


1/2 





1/2 


3 


1 




1 






3/4 


1/4 





1/4 




1/3 


2/3 





2/3 



Table 1. Parameters for Runge-Kutta time marching schemes. 



4. Fully Discrete Schemes 

Let U n is the discrete solution at t n , and At = t n+1 — t" . Then a semi-discrete scheme (|3.2p 
can be written as, 

dU n ■ 

(4.1) -^ = A, i (U") + S ili (U"), 
where, 

We describe two different time discretization schemes below. 

4.1. Explicit Schemes. We use explicit Runge-Kutta (RK) time marching schemes for the 
time-discrctizing of the two-fluid equations. For simplicity, we restrict ourselves to the second- 
and third-order accurate RK schemes (see [IB]). These methods are strong stability preserving 
(SSP). In order to advance a numerical solution from time f" to t n+1 , the SSP-RK algorithm is 
as follows: 

1. Set U(°) = U". 

2. For m = 1, k + 1, compute, 

m— 1 

1=0 

3. Set U»+ 1 = U£ +1) . 

The coefficients a m i and /3 m ; are given in Table [T] 

4.2. IMEX-RK Schemes. As discussed in section [TJ two-fluid equations contain the following 
parameters: the speed of light, mass ratio of ions to electrons, Debye length, and the Larmor 
radius. These parameters determine the time scales of the flow and may impose severe restrictions 
on the time step of explicit time marching schemes. Hence, we consider IMEX methods in this 
section. An Implicit-Explicit Runge-Kutta (IMEX-RK) scheme for (|1.1[) . is based on the implicit 
treatment of the stiff source term and an explicit treatment of the convective flux terms. This 
allows us to overcome stiffness due to the source terms. 

We will use SSP-RK schemes, as described above, with each intermediate Euler update being 
carried out by solving, 

(4.2) U™ +1 - U™.+AfAj(U m ) + AtS iJ (U m+1 ), 

for U m+1 . Usually (|4.2I) is solved using some iterative methods. However, we can exploit the 
special structure of the source term for the two-fluid equations to solve (|4.2I) exactly. We proceed 
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as follows: 

Denote U = {W 1 ,W 2 ,W 3 } with, 

W 2 = {p t v x ,p t v y ,p l v z ,PeV x ,p e v y e ,p e v z e ,E x ,E y ,E z } T 
W 3 - {E h E e ,<p} T 

We observe that (|4.2[) can be rewritten in the following three blocks, 

(4.3a) w[' Tl+1) = G!(U (m) ), 

(4.3b) W^ m+1) = G 2 (uW)li(wS m+1) )wl ra) , 

(4.3c) W^' Tl+1) = G3(U( m ))+H(W< m+1) ,w( m+1) ). 

Here Gi, G2 and G3 are the explicit parts of (|4.2j) for the variables Wi, W 2 and W3 respectively. 
The Eqns. (|4.3I) are then solved in sequential manner: 

I) Equation (|4.3al) is updated explicitly, as it involves the evaluation of the known terms from 

the previous time step. 
II) Note that the matrix _4(W^ m+1) ) in Eqn. (Obi) is, 

(4.4) 

B - - B \ ' ^ 

r g Tg r g 



,(m+l) 



S—^, — •—<. 

(m+1) 





«,(m+l) 



2, (m + 1) y ,( m+ l) 



„("> + !) 



-2-^ — --2^+ — i^, 



, Zl (m + l 



(m + 1) 



--S-4: — — ^ 



R y.(m+1) R a;,(m + 1) 

B±-i — — - B - 



„(m + l) 





^jp- 

^+ =jf- 

with f eig = —f g /\ m and i'C = X 2 r g . All the quantities in the matrix are already computed 
in step I. So, we can rewrite Eqn. (|4.3b|) as, 

(4.5) W< m+1) = (l - (At)^(W< m+1) )) G 2 (U(™)). 

which can evaluated exactly. 
Ill) The Eqn. (l43cl) is now updated for W™ +1 by evaluating H(W™ +1 , W™ +1 ). 

Remark 4.1. The IMEX scheme proposed above does not require any non-linear Newton solves 
or any global matrix inversions. It only needs explicit evaluations of the inverse of a local 9x9 
matrix in each cell making this scheme computationally inexpensive. Furthermore, there are no 
local linearizations or approximations being used in the scheme. It uses an exact solution of the 
time stepping update (|4.2I) . 

Remark 4.2. Note that the wave speeds of the system depend on the speed of light and the 
sound speeds of the electron and ion. The speed of these waves is either specified or determined 
by the flux terms of the two-fluid equations. Consequently, an explicit in time, evaluation of the 
flux terms, as in an IMEX scheme, might still lead to severe time step restrictions on account of 
these waves. 
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5. Numerical Results 

We present a set of numerical experiments to demonstrate the robustness of the proposed 
schemes. 




(a) L 1 order of convergence: L 1 -errors of the ion-density at time t 
400 and 800 cells 



2.0 are plotted for 100, 200, 



Ion-density at time t=2.0 using 1 00 cells 




2.8 - 
2.6 - 
2.4 - 




o Ref 


-Solution 




02- 


FVM-exp 




02- 


ESMinMod- 


exp 


02- 


FVM-IMEX 




02- 


ESMinMod- 


IMEX 



(b) Ion-density plots at time t = 2.0 using 100 cells 



Figure 1 . Errors of second order schemes 



5.1. Convergence Rates. As it is not possible to obtain explicit solution formulas for the 
two-fluid equations, we employ a forced solution approach to manufacture explicit solutions. 
In one space dimension, we consider the modified equation: 

d t u + d x f(u) = s(u) + K(x,t). 

with forcing term: 

K(x, t) = {0i 3 , -(2 + sin(27r(> - t))),Q, 0,2 + sin(27r(a; - t)), 0} T 
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The initial densities are pi = p e = 2.0 + sin(27ra;), with the velocities vf = = 1.0 and the 
pressures Pi = p e — 1-0. The initial magnetic field is B v = sin(27ra;) and the electric field is 
E z = — sin(27ra;). The computational domain is (0,1) with periodic boundary conditions. The 
ion-electron mass ratio is set to mj/m e = 2.0. 

It is straightforward to check that the exact solution is 

pi = Pe = 2.0 + sin(27r(a; — t)). 



In Figure l(a)| we have plotted the L 1 errors for the second-order schemes based on entropy stable 



fluxes with minmod (ES-MinMod) reconstruction for the spatial discretization and the second 
order SSP-RK scheme for time updated. For comparison, we have also plotted the results for 
the second-order FVM scheme based on a four wave HLL type solver with minmod limiter (02- 
FVM) . We observe that entropy stable schemes are significantly less diffusive than the standard 



FVM schemes. This is further verified by the solution plots in Figure 1(b) The entropy stable 
version of the IMEX scheme is also less diffusive than its FVM counterpart. However, we observe 
that rate of convergence for the IMEX scheme falls when we refine the mesh. This is on account 
of splitting errors (in each RK2 sub-step) for the IMEX schemes. 

5.2. Soliton Propagation in One Dimension. Soliton propagation in two-fluid plasmas are 
simulated in [9j [U [3j [2] . It is shown that ion-acoustic solitons can form from an initial density 
hump. In this section, we follow [SJ|3], to simulate solitons in one space dimension. 
Initially, the plasma is assumed to be stationary with ion density, 

(5.1) pi = 1.0 + exp(-25.0|ai - L/3.0|), 

and mass ratio mi/m e = 25,, on the computational domain D = (0,L) with L = 12.0. Electron 
pressure is p e — 5.0pi with an ion-electron pressure ratio of 1/100. Normalized Debye length 
is taken to be 1.0. Periodic boundary conditions are used. We consider three different Larmor 
radii: 10~ 2 , 10 -4 and 10 -6 . Numerical solutions are computed using 1500 cells. The simulations 
are carried out using an MPI parallelized version of the code, on 10 computational cores. 

The solutions are plotted for second order, spatially accurate entropy stable schemes (ESMN), 
using second (02-ESMN) and third order (03-ESMN) SSP Runge -Kutta, explicit and IMEX 
time stepping routines. We have also plotted the corresponding FVM solutions. The reference 
solutions for these simulations are computed using the 03-ESMN-IMEX scheme on 10000 mesh 
points. 

In Figure [5J we have plotted solutions corresponding to the Larmor radius of 10 -2 . This 



corresponds to the simulation performed in [5]. In Figure 2(a) we have plotted the ion-density 
profile at non-dimensional times t = 1,2,3,4 and 5. We observe that all the schemes are able 
to capture soliton waves. In particular, the speed of soliton propagation is the same for all the 



schemes. In Figure. 2(b) we have plotted the solutions at non-dimensional time t = 5.0 and 
compared them with the reference solution. We again observe that the entropy stable schemes are 
more accurate than the corresponding FVM schemes. However it is hard to distinguish between 



some schemes in Figure. 2(b) as solution lines for 02-FVM-exp, 03-FVM-exp, 02-FVM-IMEX 



and 03-FVM-IMEX coincide. Similarly, solution lines for 02-ESMN-exp, 03-ESMN-exp, 02- 



ESMN-IMEX and 03-ESMN-IMEX lie on top of each other in Figure. 2(b) To, further analyze 



the schemes in Figure. 2(c) we have zoomed in on the solution at x = 1.35. We notice that 
ESMN-IMEX schemes are slightly more diffusive than the ESMN-exp schemes. 

Compared to the solutions presented in [5], entropy stable schemes appear to be more diffusive. 
However, in [5] authors have used a fourth order Runge-Kutta update for the source updates. 
Additionally, observe that both entropy stable schemes and wave propagation method fails to 
capture the oscillation around x — 10.0 at the low resolution of 1500 cells. These oscillations are 
present in the solution only at finer resolutions. 
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Solution at time t=1 , 2, 3, 4 and 5 




- 02-FVM-exp 

- 02-ESMinMod-exp 

- 03-FVM-exp 

- 03-ESMinMod-exp " 

- 02-FVM-IMEX 

- 02-ESMinMod-IMEX 

- 03-FVM-IMEX 

- 03-ESMinMod-IMEX ' 



(a) Ion-density evolution: Ion-density at non-dimensional t = 1,2,3,4, and 5 for various schemes 

Comparison of ion-density at time t=5 



1.25 - 




Reference Solution 

- 02-FVM-exp 
-02-ESMinMod-exp 

- 03-FVM-exp 
-03-ESMinMod-exp 

- 02-FVM-IMEX 

- 02-ESMinMod-IMEX 

- 03-FVM-IMEX 

- 03-ESMinMod-IMEX 



(b) Ion-density at non-dimensional time t = 5.0 for various schemes 




(c) Ion-density at non-dimensional time t = 5.0 for various schemes: Zoomed at x = 1.35 



Figure 2. Soliton propagation using 1500 cells and Larmor radius r g = 10 
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Solution at time t=1, 2, 3, 4 and 5 





02-FVM-exp 

02-ESMinMod-exp 

03-FVM-exp 

03-ESMinMod-exp - 

02-FVM-IMEX 

02-ESMinMod-IMEX 

03-FVM-IMEX 

03-ESMinMod-IMEX - 


A 


- ^ i 


V ^ 




i 




K i 


I A 





(a) Ion-density evolution: Ion-density at non-dimensional t = 1,2,3,4, and 5 for various schemes 

Comparison of ion-density at time t=5 



1.25 - 




Reference Solution 

- 02-FVM-exp 

- 02-ESMinMod-exp 

- 03-FVM-exp 

- 03-ESMinMod-exp 

- 02-FVM-IMEX 

- 02-ESMinMod-IMEX 

- 03-FVM-IMEX 

- 03-ESMinMod-IMEX 



(b) Ion-density at non-dimensional time t = 5.0 for various schemes 



FIGURE 3. Soliton propagation using 1500 cells and Larmor radius f g = 10 



In Figures [__ and _U we have plotted the solutions for Larmor radii of 10~ 4 and 10 -6 , re- 
spectively. In Figures 3(a) and 4(a)| we have plotted the ion-density at non-dimensional times 
t = 1, 2, 3, 4 and 5. As in the previous case, we observe that all schemes capture soliton waves. 
Furthermore, from the solution plots at non-dimensional time t = 5.0, in Figures |3 (b) | and |4 (b) | 
we again note that the entropy stable schemes are less diffusive than FVM schemes. For the case 
of Larmor radius 10 -6 , we have not presented the solution for second order explicit time updates 
due to the very large simulation times, required for these schemes. 



The above figures show that the IMEX schemes are slightly more diffusive than the explicit 
schemes for the same resolution and for the same spatial discretization. A natural question that 
arises in this context is why should be IMEX schemes be used when they only differ marginally in 



16 



HARISH KUMAR AND SIDDHARTHA MISHRA 



Solution at time t=1 , 2, 3, 4 and 5 





03-FVM-exp 

03-ESMinMod-exp 

02-FVM-IMEX 

02-ESMinMod-IMEX " 

03-FVM-IMEX 

03-ESMinMod-IMEX 


A 








A - 


K A A 


I I I I I 



° 1, 



(a) Ion-density evolution: Ion-density at non-dimensional t = 1,2,3,4, and 5 for various schemes 

Comparison of ion-density at time t=5 



1.25 - 




Reference Solution 

- 03-FVM-exp 
-03-ESMinMod-exp 

- 02-FVM-IMEX 

- 02-ESMinMod-IMEX 

- 03-FVM-IMEX 

- 03-ESMinMod-IMEX 



(b) Ion-density at non-dimensional time t = 5.0 for various schemes 

Figure 4. Soliton propagation using 1500 cells and Larmor radius f g = 10~ 6 . 



Scheme 


f g = 10-' z 


fg = ID" 4 


fg = 10-" 


02-ESMinMod-exp 


100.42 


5089.67 




03-ESMinMod-exp 


152.26 


533.85 


74159.3 


02-ESMinMod-IMEX 


103.67 


106.53 


102.87 


03-ESMinMod-IMEX 


151.83 


152.3 


151.71 



Table 2. Comparison of simulation times of the numerical schemes for Larmor 
radii of 10~ 2 , 10~ 4 and 10~ 6 using 1500 cells. 
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resolution with the explicit time stepping schemes ?. The answer to this lies in the computational 
run-time. 

From the source term for the two-fluid equations (|1.5[) . we see that the Larmor radius is a 
crucial parameter in determining the strength of the source term. Reducing the Larmor radius 
leads to an increase in the strength (and hence, stiffness) of the source term. Furthermore, 
the Larmor radius does not determine the speed of the waves in the two-fluid system. Hence, 
reducing the Larmor radius is a good test for evaluating the relative advantage of IMEX schemes 
over explicit time marching schemes. 

To this end, we consider soliton propagation with different Larmor radii of 10 -2 , 10 -4 and 
10~ 6 , respectively. As the Larmor radius does not influence the wave speed, the time step for the 
IMEX schemes remains similar for the three simulations (with different Larmor radii). On the 
other hand, the increase in the strength of the source term, due to the decrease in the Larmor 
radius, implies a reduction in the time step for an explicit scheme. Therefore, we expect to see 
a difference in the computational cost between the implicit and explicit schemes. 

The simulation run-time for the three simulations (with different Larmor radii), on a mesh 
of 1500 points, with all other simulation parameters being constant, are shown in Table [5J The 
table shows that the runtime for explicit schemes increases dramatically as the Larmor radius is 
reduced. The second-order scheme is particularly affected as the stability region for RK2 is quite 
small and it requires smaller time steps. In fact, for the Larmor radius of 10~ 4 , the second-order 
(in time) explicit scheme is about 10 times slower than the third-order (in time) explicit scheme. 
As a consequence, the run-time for the second-order explicit scheme on a Larmor-radius of 10 -6 
is too large and the run was not completed. The run-time for the third-order explicit scheme 
was also very large. On the other hand, the time taken by the implicit schemes (for both second- 
and third-order time stepping) is constant with respect to the Larmor radius. This implies a 
massive speed up of the IMEX schemes (about a factor of 500) when compared to the explicit 
schemes. This example illustrates the main advantage of the IMEX schemes: their robustness 
with respect to very low Larmor radii. 



5.3. Generalized Brio-Wu Shock tube Problem. 

problem are 



The initial conditions for this Riemann 



(5.2) 



U lcft = < 



1-5 



= 1.0 

= 5 x 10" 
= 1.0 m e /mi 
= 5 x 10~ 5 



B x = 0.75 
B y = 1.0 



U 



right 



= 0.125 

= 5 x 10~ 6 
= 0.125 m e /rrii 
= 5 x 10~ 6 



E 



t/j = B z 



= 





B x = 0.75 

By = -i.o 



= v e 



E 



B z 



on the computational domain (0, 1) with, U = Uj c ^ for x < 0.5 and U 



= 


= U r jgjjj for x > 0.5. 

The ion-electron mass ratio is taken to be mi/m e = 1836. The initial conditions are non- 
dimensionalized using po = 10~ 4 . Non-dimensional Debye length is taken to be 0.01. Simulations 
are carried out using Larmor radii of 10 and 0.001. Neumann boundary conditions are used. 

The purpose of this numerical experiment is to demonstrate the behavior of the solutions of 
two-fluid equations in two different regimes: one with high Larmor-radius and another with very 
low Larmor radius, respectively. 

Numerical solutions are presented in Figure [5J In Figure 5(a)| we have plotted the numerical 
solutions based on 02-ESMinMod scheme using second order explicit and IMEX time updates. 
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Comparison of second order entropy stable explicit and IMEX schemes 




Euler 

MHD 

02-ESMinMod-exp 

02-ESMinMod-IMEX 
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(a) Generalized Brio-Wu shock tube problem: 1000 cells were used. Numerical solutions are plotted 
for f g = 10.0 



Comparison of second order entropy stable explicit and IMEX schemes 
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X 

(b) Generalized Brio-Wu shock tube problem: 50000 cells were used. Numerical solutions are plotted 
for second order schemes with f g = 0.001 



Comparison of second order entropy stable with third order explicit and IMEX time updates 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 



x 

(c) Generalized Brio-Wu shock tube problem: 50000 cells were used. Numerical solutions are plotted 
for second order schemes with f g = 0.001 



Figure 5. Generalized Brio-Wu shock tube Riemann problem 
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Solutions are computed with non-dimensional Larmor radius of 10.0, using 1000 cells. We observe 
that solution is very close to the solution of the Euler equations for each species. Note that letting 
f g — > oo, one recovers the uncoupled equations of gas dynamics for both species. Furthermore, 
both IMEX and explicit schemes produce very similar results. 

The second regime that we investigate is to set the Larmor radius to 10~ 3 . One expects to 
recover the MHD limit for vanishing Larmor radius. This limit is quite complicated to compute 
as one has to resolve the small-scale Langmuir oscillations, necessitating very fine meshes (see 
[§])• We show results obtained on a mesh of 50000 cells both for second-order and third-order 
(in time) entropy stable (explicit as well as IMEX) schemes in figures 5(b) and 5(c) respectively. 

We observe that the both explicit and IMEX solutions are converging to the MHD limit. 
However the second-order (in time) explicit scheme produces some small scale oscillations (near 
the left boundary). These small scale oscillations are not present in the results present in [5] 
as the source term in [3] is discretized using a fourth order Runge-Kutta update. On the other 
hand, the IMEX schemes resolves all the waves correctly. For the explicit schemes, small scale 
oscillations disappear when SSP-RK3 time update is used (see Figure 5(c)) and the results are 
comparable to those present in [5] in this case. 



5.4. Soliton Propagation in Two space dimensions. Two dimensional soliton simulations 
were presented in [5] . We follow [2] and simulate 2-d solitons by considering the initial ion-density 
to be 

(5.3) pi = 1.0 + 5.0cxp(-500.0((a; - L x /2.0) 2 + (y - L y /2.0) 2 ) 

on the computational domain (0, L x ) x (0, L y ) with L x = L y = 2.0. All other initial conditions 
are same as in the case of one dimensional soliton propagation in section [5~21 Neumann boundary 
conditions are used to allow the waves to exit the domain without noticeable reflections. Note 
that we consider the ion-electron mass ratio of 25 as compared to the ratio of 10, considered in 
[2J. Furthermore, we use Larmor radii of 10 -2 and 10 -4 , compared to 10 — 1 , used in [2J. We 
expect dispersive waves moving outwards, similar to the one dimensional case, considered in 
section [5~21 (Also see [5]). 



Scheme 


fg = 10-' 2 


fg = ID" 4 


03-ESMinMod-cxp 


907.2 


2661.36 


03-ESMinMod-IMEX 


921.82 


939.96 



Table 3. Comparison of simulation times of the numerical schemes for Larmor 
radii of 10~ 2 and 10~ 4 , using 200 x 200 cells. 



Numerical results are presented in Figures [6] and [Jj corresponding to the Larmor radii of 10 2 
and 10~ 4 , respectively. In Figure 
time of t = 0,0.1,0.2 and 0.3 using 03-ESMN-IMEX scheme. The wave structure observed is 
similar to the one dimensional case. In Figure |6(b)| and |7(b)| we have plotted one dimensional 
cuts of the solution at x = 1 and at non-dimensional time t = 0.5 for 03-ESMN-exp and 03- 
ESMN-IMEX schemes. As seen in the figures, the initial density hump bearks into a standing 
wave, centered at the origin, together with dispersive waves that propagate outward. We observe 
similar performances for both schemes. Furthermore, the IMEX scheme is faster than the explicit 
scheme for the low Larmor radius simulation. 



6(a) and 7(a)|we have plotted the solution at non-dimensional 
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(a) Ion-density evolution: Ion density pi at time t = 0.0, 0.1, 0.2 and 0.3. 
Ion density at non-dimensional time t=0.3 



03-ESMinMod-exp 




(b) Comparison of the schemes: Ion-density cut at x = 1 for time t = 0.3 of 
entropy stable schemes using third order explicit and IMEX time update. 

Figure 6. Soliton propagation in two dimensions on 200 x 200 mesh with r g = 10~ 2 . 

6. Conclusion 

We consider the two- fluid plasma equations and design finite difference schemes to approximate 
them. The semi-discrete version of the scheme is shown to be entropy stable. As the source terms 
in the two-fluid equations can be stiff, we propose IMEX schemes that treat the source terms 
implicitly. The novelty of our approach, in this context, is to observe that the special structure 
of the two-fluid plasma equations allows us to write the implicit (in source) time update as a 
local (in each cell) linear system of equations. This system can be solved exactly. Hence, our 
IMEX scheme does not require any Newton iterations or linearizations. 

Both the explicit and IMEX entropy stable schemes are shown to perform robustly on a set 
of numerical experiments. The entropy stable schemes are clearly more accurate than standard 
HLL type finite volume schemes. The main advantage of the IMEX schemes lie in the fact that 
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(a) Ion-density evolution: Ion density p; at time t = 0.0,0.1,0.2 and 0.3. 

Ion density at non-dimensional time t-0.3 
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(b) Comparison of the schemes: Ion-density cut at x = 1 for time t = 0.3 of 
entropy stable schemes using third order explicit and IMEX time update. 



Figure 7. Soliton propagation in two dimensions on 200 x 200 mesh with r g = 10 



they are robust (in run-time) with respect to a decrease in the Larmor radius. In particular, 
on (realistic) low Larmor radii simulations, the IMEX schemes can gain orders of magnitude in 
speedup as compared to the explicit schemes. 

We will extend the entropy stable schemes to even higher order of accuracy and couple them 
with adaptive mesh refinement to be able to simulate realistic two-dimensional examples like 
magnetic reconnection, in a forthcoming paper. 
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